Evaporation/condensation transition of the two dimensional Potts model 

in microcanonical ensemble 



Tomoaki NogawaQ and Nobuyasu Ito 
Department of Applied Physics, The University of Tokyo, Hongo, Bunkyo-ku, Tokyo 113-8656, Japan 

Hiroshi Watanabe 
Institute for Solid State Physics, The University of Tokyo, 
Kashiwanoha 5-1-5, Kashiwa, Chiba 277-8581, Japan 

Evaporation/condensation transition of the Potts model on square lattice is numerically investi- 
gated by the Wang-Landau sampling method. Intrinsically system size dependent discrete transition 
between supersaturation state and phase-separation state is observed in the microcanonical ensem- 
ble by changing constrained internal energy. We calculate the microcanonical temperature, as a 
derivative of microcanonical entropy, and condensation ratio, and perform a finite size scaling of 
them to indicate clear tendency of numerical data to converge to the infinite size limit predicted by 
phenomenological theory for the isotherm lattice gas model. 



I. INTRODUCTION 

First order transition as a state transformation of sub- 
stances is observed in our everyday life, such as evapo- 
ration/condensation and melting/freezing. The thermo- 
dynamic mechanism of such transition is very clear; a 
phase transition occurs when the state of the lowest free 
energy alternates from one thermodynamic state to an- 
other by changing a certain environmental parameter. 
But dynamics of the transition is rather complicated and 
our knowledge cannot be said to be sufficient. Since 
the initial and final state is completely different unlike 
a second order transition, nucleation and large scale do- 
main growth, i.e., invasion of the metastable state by the 
most stable state, is observed in the vicinity of the tran- 
sition point, which includes various mechanisms depend- 
ing on the spatio-temporal scale [1 - 3] . This is essentially 
nonequilibrium phenomena and the dynamics is difficult 
to understand by relating it to the well-understood equi- 
librium state near the transition point. There is a droplet 
formation phenomena, however, that can be discussed in 
an equilibrium framework as mentioned in the following. 
It can be a good starting point to understand first order 
transition dynamics. 

About first order transition, we usually imagine discon- 
tinuous transformation with hysteresis controlled by in- 
tensive variable such as temperature, pressure and mag- 
netic field. But once we constrain one conjugate exten- 
sive variable, such as internal energy, particle density and 
magnetization, and take it as a control variable, a coexist- 
ing phase is inserted in the phase diagram between two 
homogeneous phases and the transitions becomes con- 
tinuous; the volume fraction of one phase continuously 
changes from zero to unity between the two edges of 
the coexisting phase. It is pointed out, however, that 
a droplet, which is a condensed domain of the minority 
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state in phase separation, is not stable unless its volume 
is larger than a certain threshold and therefore there is a 
discontinuous droplet condensation transition [4, 5] . This 
is basically due to that the surface free energy of the 
droplet is not negligible in a finite size system. Therefore 
the threshold volume of <i-dimensional droplet depends 
on the system size, proportional to L d where L 

is a linear dimension of the system [||. Although this 
threshold can be neglected in comparison with the whole 
system volume L d in an infinite size system, it actually 
diverges with L. Consequently, the infinite size limit of 
the droplet condensation transition is well-defined if we 
observe the phenomena with proper size scale. It should 
be noted that the equilibrium droplet condensation in a 
large size system always contain single droplet in contrast 
with nonequilibrium transition where multiple droplets 
appear. The supersaturation state, where condensation 
is avoided due to small volume of minority state, is re- 
lated to the metastable state under fixed intensive pa- 
rameter condition. 

Biskup et al. Q and Binder Q made quantitative anal- 
ysis of the equilibrium droplet condensation transition 
which is supported by a number of numerical studies of 
Lenard-Jones particles 0-0 and lattice gas fl(il - [l3} where 
particle density is constrained at given temperature. In 
this paper we consider a simpler situation where rele- 
vant extensive variable is only internal energy (magneti- 
zation is not taken account of) and treat the equilibrium 
state with constrained energy, i.e., microcanonical en- 
semble [HI, EH ■ By large scale numerical simulations of 
the two dimensional Potts model with the Wang-Landau 
sampling, we try to determine the large size limit of the 
droplet condensation transition, which has been rather 
difficult in small size systems flOf . 



II. MODEL 

We investigate the ferromagnetic g-state Potts model 
flfjj on L x L square lattice with a periodic boundary con- 




FIG. 1. Typical spin configuration of the system with q = 8 
(top) and q — 21 (bottom) for L = 256. The different spin 
states are denoted by different gray scales. The energy is 
E/N = 0.3750,0.4500,0.5250 and 0.7125 (top) and E/N = 
0.1500, 0.2250, 0.2625 and 0.6750 (bottom) from the left to the 
right. The black region is occupied by spins of the most major 
state and dappled gray region is assemblage of tiny domains 
with various states. 

dition. The interaction energy for the nearest neighbor 
pairs is written as 

e= £ (i-<W, (i) 

(i.j)en.n. 

where 8 a p means Kronecker's delta and the spin variable 
Cj takes integer value 1, 2, • • ■ , q. This model in canonical 
ensemble exhibits paramagnetic-- ferromagnetic transition 
at P = p c = ln(l + l/y/q) where f3 = l/k B T is the in- 
verse temperature. The transition is of second order for 
q < 4 and of first order for q > 4. When q is not suf- 
ficiently larger than 4, the correlation length of fluctua- 
tion is considerably large around the transition point and 
critical like fluctuation is observed to some extent. It is 
better to choose large q to investigate the pure nature of 
a first order transition. But the amount of computation 
for equilibration becomes larger with increasing q. 
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FIG. 2. Schematic diagram of dividing energy region to 8 
threads and spin configuration exchange. 
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FIG. 3. (color online) The standard deviation of microcanon- 
ical temperature among samples with different realization of 
random number. The data without exchange and 4 and 16 
exchanges in every Monte-Carlo steps per spin are plotted to- 
gether. The simulation condition is q = 8, L = 512, 16 sam- 
ples and the modification constant is down to In g(E) = 2 -20 . 
The energy region, 0.36 < E < 0.56, is distributed to 64 
threads and each threads handle the energy points about 
666 x 2. Smoothing is done by averaging for 128 energy 
points after calculating deviation. For comparison, the dif- 
ference between the transition temperature and the spinodal 
temperature, a characteristic scale of our interest, is about 
0.001 for this size as shown in Fig. [4ja). 



III. METHOD 

We p erform the Wang-Landau sampling simulations 
fiTL 1 1 8j] - which yields the density(number) of states g(E), 
as a result of learning process to realize a flat energy 
histogram. The density of states enables us to calcu- 
late the Hclmholtz's free energy in canonical ensemble as 
F(/3) = ^ _1 ln [Ei?. 1 ?^) 6 "^]) and its derivatives, i.e., 
mean energy and specific heat. The energy with max- 
imum or minimum realization probability is given as a 
solution of d[g{E)e-^ E ]/dE = 0, i.e., din g(E)/dE = /3. 

On the other hand, in microcanonical ensemble for 
given internal energy, fundamental thermodynamic func- 
tion is microcanonical entropy S(E) = fcslng(i?) and 
(inverse) microcanonical temperature is a quantity to be 
observed, which is defined as a response to energy per- 
turbation as 



1 dS(E) _ d_ 
dE ~ dE 



In «?(£). 



(2) 



This is equivalent to the extremal condition for the free 
energy in canonical ensemble. The independence of /3 
in microcanonical ensemble is not exactly related to the 
/3 dependence of expectation value of E in canonical en- 
semble except in the thermodynamic limit but is exactly 
related to the /3-dependence of the extremal value of E 
even in the finite size system. 

We perform parallel computation to treat large size 
system as done in Ref. [IH . The energy region is divided 
into a number of parts with constant width and each 
part is associated to a different thread. The spin flip tri- 
als which make the energy of the system go out of the 
given range is always rejected. Since the time needed to 
diffuse the energy range AE by random walk is propor- 
tional to AE 2 , the time needed to obtain flat histogram 
is inversely proportional to the square of the number of 
threads. Although this method drastically reduces the 
total Monte-Carlo steps, we have to care about the possi- 
bility that a flat histogram is established in Monte-Carlo 
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steps shorter than the relaxation time of the system. This 
occurs when the energy region is divided into too small 
parts. Another problem of the division of energy region 
is that it possibly causes the segmentation of the phase 
space, i.e., there are spin configurations which cannot be 
visited depending on the initial condition. 

In order to enhance the relaxation and guarantees er- 
godicity, we make overlapping energy region for neighbor- 
ing threads as illustrated in Fig. [21 where the exchange of 
spin configuration is allowed satisfying a detailed balance 
condition, 

W({X,Y} -> {F, A}) = 9i {E{X))g 3 (E(Y)) 
W({Y, X} {X, Y}) 9l (E(Y)) g 3 (E(X)) ' 

where X and Y are indices of microscopic states, 
in the same spirit with the replica exchange method 
fl9j . Here, gi is a density of states calculated by 
i-th thread and W({X,Y} -> {Y,X}) means the 
transition probability from a compound state: X 
for the i-th thread and Y for the j-th thread, to 
its exchanged state. Practically the exchange is ac- 
cepted with a probability, W({X,Y} -> {Y 7 X}) = 
min(l, 9l (E(Y)) g 3 (E{X)) / 9i (E(X)) g 3 (E(Y))), simi- 
larly with the Metropolis method. 

In order to check the efficiency of the replica exchange 
method, the sample to sample deviations of temperature, 

5/3 = \J ft 2 — |/3| 2 , where 7TT means average over samples, 
are shown in Fig. [3] The conditions of simulations are de- 
scribed in its caption. The peaks arc observed in the sys- 
tem without exchange. The amplitudes of the peaks be- 
comes smaller as the exchange frequency increases. The 
peaks exist at the seam points, which are the edges of 
the divided region. Except around the seam points, the 
reduction of fluctuations are moderate. This is because 
the present system does not exhibit large fluctuation for 
microcanonical state. The exchange method can work 
more effectively for the system with more complicated 
landscape in phase space. 



IV. RESULTS 

We perform simulation with q = 8 and 21 for system 
size L = 32 — 1024 by using 4 — 64 threads. The density of 
states is calculated not for all energy region but restricted 
region of interest. The numerical data shown below are 
average over 2 — 8 samples with different realizations of 
random number. In the Wang-Landau sampling we de- 
crease the modification constant for \n g(E) step by step, 
l,2-\2- 2 ,---,2- 30 (except for q = 8 with L = 1024 
and q = 21 with L = 512, where modification constant is 
decreased down to 2~ 25 ) with achieving 90 % flatness of 
energy histogram for all threads in every step. Typical 
number of total MCSs is 10 6 — 10 7 per spin in the final 
step of simulations. The replica exchange is attempted 
in every 16 MCSs per spin. 
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FIG. 4. (color online) Inverse temperature as a function of 
internal energy for (a)q = 8 and (b)q — 21. Smoothing is 
performed by averaging over the range y/~N /4 of E. The hor- 
izontal line indicates fi — ln(l + 1/^/q). The vertical lines 
indicates E/N = eZ (left) and et (right). The arrows indi- 
cates the direction that L becomes larger. 



A. temperature vs internal energy characteristics 

Figure|4]shows the inverse microcanonical temperature 
(3(E) = g(E + 1) — g{E) for q = 8 and 21 for various sys- 
tem size. While f3(E) should be a monotonically decreas- 
ing function of E to make the free energy a convex func- 
tion of E, it is not for the finite size system. In the region 
with positive derivative, i.e., negative specific heat, phase 
coexisting is observed as shown in Fig. [T] This state cor- 
responds to the free energy maximum in the canonical 
ensemble. A dip and a peak exist inside the coexisting re- 
gion in the thermodynamic limit, e~ < E/N < e+, where 
e~ ~ 0.403 and e+ ~ 0.888 for q = 8 and e~ ~ 0.171 
and e+ ~ 1.392 for q = 21. We note the bottom/top 
position of the dip/peak as (E,f3) = (E~ pi (L) , (3~ pi (L)) 

and (£^;(L), (3^ pi (L)), respectively. These points are re- 
garded as the equilibrium spinodal points, i.e., saddle- 
node bifurcation points, where the second free energy 
minimum in canonical distribution annihilates together 
with one free energy maximum. As increasing system 
size, E^pJN approaches ef and /3 approaches f3 c for all 
E with e" < E/N < e+ 
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FIG. 5. (color online) Finite size scaling results of the excess 
temperature. (a)g = 8 and (b)q = 21. The arrows indicates 
the direction that L becomes larger. The thick gray curve is a 
guide for eyes. The vertical line in the bottom panel indicates 
E-E~ = 1.507L 4/3 as well as in Fig. E 
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FIG. 6. (color online) Probability distribution function of the 
energy density of subsystems for q — 21, L = 512. Each 
curves corresponds to E/jV=0.12-0.52 (0.02 step). 

approaches a large size limit; |/3 — f3 c \ linearly increases 
as (k B (3 2 /C ± )\E — Ef \ where is a specific heat, 
-k B P 2 [d 2 \ng(E)/dE 2 }~ 1 at E = Ef in the evaporation 
phase and discontinuously drops down when entering the 
condensation phase. Note that this scaling does not aim 
the collapse of data to a universal curve as in standard 
finite size scaling for second order transitions, but it is 
intended to show the conversion to the large size limit 
by blowing up the transition region, which vanishes in 
macroscopic scale. (Scaling of finite size rounding may 
be a challenging problem.) On the other hand, discontin- 
uous behavior is rarely observed in the system with q = 8 
even for L = 1024. This is because correlation length is 
rather large and the used sample size is effectively much 
smaller than that for q = 21. 



For q = 21, a plateau, where (3 almost equals /3 C , is 
obviously observed in the middle of coexisting region. 
In this plateau region, the two coexisting phases form 
strips parallel to the horizontal or longitudinal axis, while 
minority phase forms a droplet in the side regions (see 
Fig. [TJ. With increasing system size, it is observed 
that the plateau region becomes wider and the step-like 
change of (3(E) on the edge becomes sharper as a conse- 
quence of droplet-strip (slab) transition [10| . 

As increasing system size, the two spinodal points, 
(E^(L),^{L)), approach (£±,/3 c ) = (JVe± re- 
spectively. The deviation \^ pi (L)-(3 c \ and \E± pi (L)-Ef\ 
decreases as a power function of L [2(| • To show it clearly, 
we perform finite size scaling in Fig. [5] with expecting a 
formula 

\m-^=L-^ F {^^), ( 3) 

with a scaling function F(-), that for isotherm lattice 
gas model is derived in Ref. Q, where /3 and E are re- 
placed with magnetic field and magnetization, respec- 
tively For q — 21, it is observed that the scaled curve 



B. volume fraction of a droplet 

To observe the nature of droplet formation more di- 
rectly, we estimate the ratio of condensation volume to 
the whole system volume. To this end, we divide a sample 
into L square-shaped subsystems with size 
is approximated by the closest integer if necessary) . This 
resolution is fine enough to capture the shape of the criti- 
cal droplet with size oiO{L d ^ d+v >). We calculate energy 
per spin el for each square, whose distribution function 
P(el) is shown in Fig. O Hereafter we only show the re- 
sults for q = 21. Bimodal distribution is observed in the 
condensed regime, E~ pi (L) < E < E+ pi (L). The width 
of peaks becomes narrower as y/sZL/L ~ L -1 / 2 and the 
integral of bridge component between the two peaks cor- 
responding to the perimeter of droplet becomes smaller 
with (2nL x y/Z)/N = L~ 1/2 . The positions of the 
peaks hardly change with E but only heights change for 
E~ pi (L) < E < E^ pi (L), which means that the change of 
state in this regime can be described only by the change 
of mixing ratio of the two phases. We determine a sub- 
system with energy below/above £ m = (e+ + e~)/2 = 
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FIG. 7. (color online) Internal energy dependence of the vol- 
ume fraction of the high energy domain. The insets are the 
blowups of the edges of the coexisting phase, where the arrows 
indicates the direction that L becomes larger. 



(1 — l/y/q) [3 belongs to the ferro/para domain. 

The volume fraction of the para phase, P(el > s m ), is 
plotted in Fig. [7] In the thermodynamic limit, P(el > 
Em) equals zero for E < E~, unity for E > E+ and 
linearly increases in between as (E — E~ ) / (£"+ — E~ ) . 
Finite size deviation is observed on the edge of the coex- 
isting region as shown in the insets of Fig. [71 The fraction 
of condensed phase for finite size system is smaller than 
that for L = oo. The normalized condensation ratio 
P(e L > e m ) 



A = 



< 1 



(E — Ec)/ (Ec — E L , 
is expected to be a function of dimcnsionlcss variable, 

(E-E-)^/ d 



(4) 



A = 



with 



L d 

g £ (£ ± -£_)^^ 



(5) 
(6) 



for L — > oo. Here %• is an interface free energy per 
volume of an optimally shaped large Wulff droplet |2lj . 
The derivation of this formulae is described in Appendix 
IA"1 The behavior in the thermodynamic limit is exactly 
known as A = for A < A c = (1/2) (3/2) 3 / 2 and 
l/4\/A(l - A) = A for A > A c . While C" is estimated 
as 5.41, tw is an only unknown quantity needed to cal- 
culate A. By only assuming tw = 0.40 (then a = 0.44), 
however, our numerical result shows good agreement with 
the theory as shown in Fig. [5] While A decreases to zero 
for A < A c with increasing L, very little finite size de- 
pendence is observed A > A c . The gap of A at A c also 
coincides with the predicted value 2/3. 



V. CONCLUSION 

We have investigated the condensation/evaporation 
transition of the Potts model in microcanonical ensem- 
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FIG. 8. (color online) Finite size scaling result of the volume 
fraction of condensed droplet. The data for E/N < 0.40 are 
used. The solid gray curve indicates theoretical prediction, 
a\E - E-^ i/2 /L 2 = 1/4 V / A(1 - A) for A > 2/3, where we set 
a = 0.44. The vertical line indicates E — E~ = 1.507L 4/3 as 
well as in Fig. [5] 



blc. Interesting property of this transition is that it be- 
comes impossible to be found in the true thermodynamic 
limit; the width of supersaturation regime disappears, 
but discontinuity of scaled quantities becomes clear with 
increasing system size. The present numerical results 
with q = 21 show good agreement with the theoretical 
prediction of the system size dependence both for tem- 
perature [j| and condensation ratio [H . Both of the two 
mean that droplet with size smaller than 0(L d /( d+1 )) is 
unstable. 

For finite size system in microcanonical ensemble, we 
observe negative (inverse) specific heat in the coexisting 
region and ovcr/undcrhang of temperature, that corre- 
sponds to the thermodynamic spinodal point in canonical 
ensemble. The scaling behavior |/3 sp i — f3 c \ oc L~ d /( d+1 ^ 
suggests an existence of a diverging length scale, R B {P) oc 
\/3~ /3 c \~( d+1 ^ d . This means that a supersaturation state 
at given temperature /3 becomes unstable at a length 
scale above R s . Although this length is related to the 
equilibrium spinodal point, it is not clear whether it also 
has some meanings in noncquilibrium dynamics, which 
may be an interesting open problem. 

Last of all, let us note on the difficulty of the simula- 
tion of first order transitions. We suffered from slow re- 
laxation for large system size in the present study as well 
as discussed in Ref. [l(| ■ It is considered that the Wang- 
Landau sampling is quite efficient for first order transi- 
tions with discontinuity of 0(L d ) as well as some other 
extended ensemble methods, such as the multicanoni- 
cal method [22| because it provides additional proba- 
bility weight on coexisting states to bridge the divide 
between two distinct homogeneous states, such as para 
and ferro phases. However there are still discontinuous 
transitions in the coexisting phase, that is a evapora- 
tion/condensation transition with 0(L d /( d+1 )) discon- 
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tinuity and droplet/slab transition, [H|- This be- 
havior is general for first order transitions. Since the 
problem is that the two states are energetically degener- 
ate in these transitions, one fundamental solution to this 
problem may be to employ another argument which cor- 
responds to the shape of domains for the joint density of 
state [3 [H[ . Although it requires considerable amount 
of computation, massive parallel computation makes it 
feasible. 



Appendix A: finite size effect on droplet 
condensation ratio 



T W V, 



(d-l)/d 



A(l - A) 2 + A (d+1 H (A2) 



with A 



- e-) 



(d-i)/d 



(6E) 



(d+l)/d 



2C tw 



L d 



(A3) 



The first term indicates the fluctuation without phase 
coexistence characterized by specific heat C~ and the 
second term indicates the surface free energy of a large 
droplet. 

The nature of the fluctuation depends on the dimen- 
sionless parameter A which is related to the excess en- 
ergy. For given A, free energy F is minimized at 



Here we derive the condensation rate of a droplet for a 
energy driven phase transition by translating the result 
for a magnetization driven transition [H, HH . We consider 
canonical ensemble at the bistable point j3 = j3 c and its 
energy distribution function around a peak at E~ . (The 
case around the another peak Ej is derived in the same 
way. ) 

We note the excess energy in fluctuation beyond the 
peak as 8E = E — E~ and divide into two parts, 
SE = 5Eb + SEd, where SEb is due to small bubble ex- 
citation and SEd is due to a large droplet of disordered 
phase. If introducing condensation ratio A as SEd = XSE, 
the rest is given by SEj, = (I — X)SE. The volume of the 
droplet, Vd = AVl, can be smaller than that in the ther- 
modynamic limit, Vl = SE/(s^ — e~), for finite size sys- 
tem. By using these quantities, the energy distribution 
function P(f3 c ;E) is proportional to e~^ cF , where 



Pc{SE b f 
2L d C- 



r w V ri 



(d-l)/d 



(AI) 



A = for A < A c = (l/2)(3/2) 3/2 , 



(A4) 



which means supersaturation regime, while A for mini- 
mum F is given by a solution of 



1/4VA(1-A)=A for A>A C 



(A5) 



which means condensation regime. The two states with 
A = and A = A c = 2/3 is equivalently stable at A c . 

This discontinuous change at A c is directly observed 
as a internal energy driven transition in microcanonical 
ensemble. 
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